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Abstract 

A lattice Boltzmann (LB) scheme is described which recovers the 
equations developed by Qian-Sheng for the hydrodynamics of a ne- 
matic liquid crystal with a tensor order parameter. The standard 
mesoscopic LB scalar density is generalised to a tensor quantity and 
the macroscopic momentum, density and tensor order parameter are 
recovered from appropriate moments of this mesoscopic density. A 
single lattice Boltzmann equation is used with a direction dependent 
BGK collision term, with additional forcing terms to recover the anti- 
symmetric terms in the stress tensor. A Chapman Enskog analysis is 
presented which demonstrates that the Qian-Sheng scheme is recov- 
ered, provided a lattice with sixth order isotropy is used. The method 
is validated against analytical results for a number of cases including 
flow alignment of the order tensor and the Miesowicz viscosities in 
the presence of an aligning magnetic field. The algorithm accurately 
recovers the predicted changes in the order parameter in the presence 
of aligning flow, and magnetic, fields. 
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1 Introduction 



The lattice Boltzmann (LB) method has been extensively studied as a meso- 
scopic method of simulating isotropic fluids {eg |l|, The strengths 
of the method lie in modelling flow in complex geometries {eg 0) or in 
multi-component flow {eg @, 0). Recently, a number of LB schemes have 
been developed to represent the flow of anisotropic fluids such as liquid crys- 
tals 1,1. 

Materials which exhibit liquid crystal phases have anisometric molecules 
|llO| , p]T[] and nematic liquid crystals are of particular interest because of their 
application in display devices. There is also increasing interest in modelling 
liquid crystal colloids in which colloidal particles are embedded in a nematic 
phase The colloidal particles interact through the distortions and defects 
which they generate in the nematic elastic field. The forces between the 
particles and the dynamics of their motion is accessible to experiment (c/ 
Poulin et al |]T^). As a consequence of these additional colloidal interactions 
the particles rearrange themselves into new structures and the LB method 
provides a particularly effective way of dealing with the complex boundary 
conditions in such problems. The method can recover both the nematostatics 
and nematodynamics of these phases and hence the dynamics of the colloidal 
phase can be captured and this has been achieved with an earlier LB method 
developed by the authors [§,0]. Additionally, the coHoidal particles may be 
an isotropic fluid and work is currently in progress to extend the two phase 
LB algorithms which have been developed for isotropic fluids to describe a 
mixture in which one of the phases is a nematic liquid crystal. This current 
work is based on the algorithms presented in this paper extended to embody 
the isotropic-nematic interface following the macroscopic description of Rey 



15, 16]. The hydrodynamics associated with pair annihilation of line defects 



is also a matter of current interest and results have been reported using 



both lattice Boltzmann |17] and conventional solvers ||T8[ of the equations for 



nematodynamics in the presence of a variable order parameters. 

The orientational ordering of the molecules in the nematic phase is char- 
acterised by a director field, nQ(x, t), a unit vector which essentially defines 
the 'average orientation' of the molecules. However, the nematic ordering is 
more fully characterised by a traceless and symmetric order tensor, Qaf3- In 
this work, for simplicity, we assume that the director is confined to a two 
dimensional plane, and hence that Q^ii may be written in the form 
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Qap{x,t) = S{2nanp - 5ap) (1) 

The principal eigenvector of Qa/j is the director and the principal eigenvalue, 
^(x, t), is the scalar order parameter. 

The scheme proposed by Care et al to recover the Ericksen, Leslie and 
Parodi (ELP) equations ||10| used two coupled LB equations. These equations 



govern the continuum fluid dynamics of an incompressible nematic with an 
order parameter, S*, which is both position and time independent. In the 
scheme P , one of the LB equations carried the momentum and the second 
carried a vector density corresponding to the director fleld. The Denniston 
et al 1^ scheme is also based on two coupled lattice Boltzmann schemes one 
of which carries a momentum density and the second of which carries an 
tensor density from which the macroscopic order tensor can be recovered. 



This latter scheme recovers the Beris-Edwards equations [jT9[ for the flow of 
a nematic liquid crystal with variable order parameter. 

In the work presented in this paper, a third scheme is presented which 
recovers the Qian-Sheng [^| equations for the flow of a nematic liquid crystal 



with a variable scalar order parameter. The scheme is based on a single LB 
equation which governs the evolution of a tensor density and from which both 
the macroscopic order and momentum evolution equations are recovered. It 
is assumed that the dynamics generated by the Qian-Sheng equations are 



essentially equivalent to those generated by the Beris-Edwards equations |[T9|] 
although this equivalence has not been demonstrated explicitly. Recent work 
by Sonnet, Maffettone and Virga provides the basis upon which the 



variety of schemes with a variable order parameter may be compared. The 
Qian-Sheng equations, like the Beris-Edwards equations, reduces to the ELP 
formalism in the limit that the order parameter becomes independent of 
time and position. There are a number of important differences between the 
scheme proposed here and that of Denniston. In the scheme described here, 
the target equations are those of Qian-Sheng rather than Beris-Edwards, 
there is a single Boltzmann equation, the equilibrium distribution function 
is isotropic (as is expected on physical grounds), and the scheme is able to 
recover the full tensorial coupling of the order tensor to the velocity gradient 
tensor, from a single lattice distribution function. 

The target macroscopic equations of the Qian-Sheng scheme are sum- 
marised in Section |^, the new scheme is described in Section ^ and a Chap- 
man Enskog analysis of the scheme is presented in Section ^. A number of 
analytical results are developed from the Qian-Sheng equations in Section |^, 



3 



which provide the basis for vahdation of the LB algorithm in Section ^. The 
conclusions are presented in Section ^ 



2 Qian-Sheng formalism 



In this section we summarise the target macroscopic equations for the LB 
method. The two governing equations of the Qian scheme are the momentum 
evolution equation 

pDtup = df,{-P6a(s + <yip + <yi^ + <y'^p) (2) 
and the order tensor evolution equation 

JQaP = ha/S + h'^fS ~ ^^aP ~ ^aP^\ (3) 

It is shown by Qian that in the limit of constant order parameter, the solu- 
tions of these equations is identical to those obtained from the ELP equations. 
Throughout this work we use the repeated index notation for summations 
over Cartesian indices. In the above equations Dt = dt + u^dfj, is the con- 
vective derivative, P is the pressure, o"^^ is the distortion stress tensor given 
by 

and a^^p is the stress tensor associated with an externally applied field. In 
the current work we only consider an externally applied magnetic field for 
which the stress tensor is given by Landau and Lifshitz |^ 

^i/3 = -^(^«^/3-^^'M (5) 

The viscous stress tensor cr^^ is given by 

^Ai2iVa/3 - fJ'lQaf.Nf.f^ + ^llQp^N^oc (6) 

where Aap = \{daUp + dpUa) is the symmetric velocity gradient tensor. The 
elastic molecular field is given by 
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where, in the current work, the Landau-deGennes free energy is assumed to 
be of the form 

(3FQ,uQurQr^ + MQ%Y (8) 

where Qai3,'y = d^iQap)- Using this form of the free energy and the definition 
of the molecular field, equation (^, we find 

The quantities A and Aq, in equation (|]) are Lagrange multipliers which im- 
pose the constraints on the elastic molecular field which arise because the 
order tensor, Qap-, is symmetric and traceless. For a three dimensional sys- 
tem they have the values 

A = \{hf_ni — 1/^2^/1/^) Aq, = \sa^uhfj,u (10) 

The term A^^ in the first of equations ([T0|) does not appear in the Qian-Sheng 
equations, but is necessary in order to correct for the slight compressibility 
of the LB fluid. In the presence of an external magnetic field the free energy 
is assumed to be augmented by a term of the form 

FH = -\{x\\-XA.)Qo.pH^Hp (11) 

which gives an additional term in the molecular field of the form 

h^p = \xaH^Hp (12) 

where Xa = X\\ ~ Xi. is the anisotropy in the susceptibility. The viscous 
molecular field, /i^^, is given by 

= ^ At2-4a/3 + /iliVa/3 (13) 

where Najj is the co-rotational derivative defined by 

Na/3 = dtQa(3 + U^d^Qap — £afiu^J-'fiQi'l3 — £ I3^iu'^ tiQ ua (14) 

where the vorticity, uj_= \{V x u). 
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In order to develop a lattice Boltzmann scheme we re-arrange equation 
(P) by substituting from equation ([T^ ) for the co-rotational derivative, A'q,^, 
and obtain the form 

//il Z 

(15) 

Equation (^) for the order tensor evolution can be cast in the form 

Z/ii /il 

^(^a/jA + £:a/3^A^) (16) 

where we have set the inertial density of the fluid, J, to zero. Equations ( P^ ) 
and (0) are the expressions which are recovered by the LB scheme proposed 
below. After each time step the macroscopic density, momentum and the 
order tensor are recovered from the LB tensor density. The gradients of the 
order tensor are then used to construct the molecular field through equation 
d^) and this modifies the dynamics through appropriate forcing terms in the 
LB equation. 



3 The lattice Boltzmann algorithm 

In order to recover both the momentum evolution equation (|^) and the order 
tensor evolution equation (0) within a lattice Boltzmann scheme, the scalar 
density of a standard lattice Boltzmann scheme, fi{r,t), is replaced by a 
tensor density giajsi'Lit) where i is the normal velocity index and a and /3 
label either a two, or three, dimensional Cartesian basis. One can think 
of the tensor density as carrying information about the ordering of that 
population of the fluid 'element' on the velocity link i, associated with a 
particular position and time. Hence the densities of a standard LB scheme 
have been generalised to carry information about the order associated with 
the fluid in addition to the density and momentum. 

The density giaisiL^t) is assumed to have the following moments 

P = ^Qif^f^ (17) 
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pUa = J^Ciafi'iMM (18) 

i 

pSaf3 = "^giafS (19) 
i 

where Sajj is an order matrix with unit trace which is related to Qa/B by 

Sap = ^(<5q/3 + Sap) (20) 

where d is the dimensionahty of the space which the Cartesian indices a and 
(3 span. The LB algorithm governing the evolution of Qiap is taken to be of 
the form 

9iap{r + 6ci,t + d) = giap{r, t) + ^ 9jP\i:, t)M,o,pj^,^ + 0^^^ + Xiap (21) 

3 

where (piap and Xiap are the forcing terms for the momentum and order 
respectively and there is an implicit summation over repeated Cartesian in- 
dices. The non-equilibrium distribution function is given by 

(neq) _ (0) /^iqn 

,(0) 
ia/3' 



where the equilibrium distribution function, gl'^g, is taken to be of the form 



= S„,ff>. (23) 

The distribution function, f^^\ is assumed to be second order in the velocity 
Ua and is determined in the usual way, by the requirements 

p=j:f''' pua=Y.^.jr (24) 

i i 

and Galilean invariance in the form 

CiaCipgf^p = Clp6al3 + pUaUf) (25) 

i 

Hence, we have 

Yl & = pSaP ^^"^S = pSal3Ua (26) 



In Section ^, it is shown that in order to recover the required tensor 
coupling of the order tensor to the velocity gradient tensor it is necessary to 
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have sixth order isotropy of the velocity tensors. In this work, the required 
isotropy is achieved by using a D2Q13 three speed hexagonal lattice with 
velocity vectors, Cp, given by 



Co = {0,0} 

ci = c{±l,0},c{±l/2,±v/3/2} 
C2 = c{0,±v/3},c{±3/2,±v/3/2} 



(27) 



the subscripts 0, 1 and 2 being associated with particles with velocity 0, c 
and ^/3c respectively. Imposing the conditions (|2^) and (^Sf ) leads to the 
result 

(28) 



(0) _ c + 



, 1 1 



6, 



al3 



where, in D2Q13, to = 11/25, ti = 9/100 and = 1/300 and the velocity of 
sound Cs = (3/10)^/^. The collision operator is taken to be of the form 



OafiOfSu I 1 



(29) 



'0 'J "^s 'J 

where the direction dependent relaxation parameter is defined to be 

ro(l + r]QCi^Ciu5^y + r}2Ci^,CiyQ (30) 



Ti 



This collision operator may be seen more transparently after performing the 
contraction with the non-equilibrium distribution function as in equation 



(neq) 



(neq) 



where 



(5 m,- 



'ia/3 



^3 



(31) 



is a correction to conserve mass. The term 

5piaP -- 



tiCis -v ^ (neq) 



(32) 
(33) 

(34) 



includes a term in I/tj which is a correction to conserve momentum and a 



term associated with a factor — 2 which is explained at the end of Section O 
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The collision operator is essentially an LBGK collision operator (e(7 0) with 
a relaxation time which is direction dependent. When written in the matrix 
form, the collision operator is seen to lie between a conventional LBGK 



operator and a linearised lattice Boltzmann scheme {eg [^) although the 
usual circulant properties of the matrix associated with an isotropic fluid 
are destroyed by the direction dependent scattering in the method presented 
here. 

The momentum forcing term is given by 

(pia/S = USa(3Cifj,duFy^ (35) 

where 

1 / fi2{,ha/3 — £a/3ti^iJ,) 



^2 



and the angular forcing term is given by 

TT- [P^Kp - ^Kp + 2(5„^A + 2£:„^^A^] (37) 

4 Analysis of the algorithm 
4.1 Chapman Enskog expansion 

In this section the key results of the Chapman Enskog analysis are given 
which demonstrate how the scheme described in Section ^ can be recovered 
from the LB algorithm described in Section ^ We follow a standard Chap- 
man Enskog analysis ( 0) and expand the density and the time derivatives 
in the form 

oo 

g^ap = Y^^^gtJp (38) 

n=0 
oo 

dt = (39) 

n=0 
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We assume for the purposes of the following analysis that the forcing terms, 
(f)ia/3 and Xia/3, cau both be introduced at order 0(5^) since both include 
gradient terms in both the velocity and director field (c/ 0). The choice 
is supported by the agreement of the Chapman Enskog analysis with the 
measured simulation data; however the assumption relies upon assumptions 
implicit in LB hydrodynamics in general. Accordingly, it requires more care- 
ful analysis which will be undertaken in later work. 

The Chapman Enskog expansion gives to 0{6) 

9t,o9l% + Cif,d^gl% = 9j]}^^iafSj,.u (40) 



2\ 



and to 0{6 

X] gfju^iaf^jf^u + (Piap + Xia/B (41) 
j 

We now take moments of equations (|40|) and (|4l|) in order to recover the 
macroscopic equations to which the algorithm is equivalent. If we sum the 
two equations over the index i, take the trace of the distribution density and 
sum the resulting equations, we recover 

{dt,o + dt,i)p + d^{pu^) = (42) 

which is the continuity equation to second order. If we multiply equations 
(|40|) and ( ^8]) by Ci^, sum over the index i, take the trace of the distribution 
density and sum the resulting equations we recover 

{dt,o + d,,i){pu^) + 9,(ng) + n« + ifii;)) = c%F,^ (43) 

where 

nS = E^-^^/^4t (44) 

i 

and 

= Yl djiu^iaCipMi^^j^^ (45) 

Further progress can be made by noting that, to 0{u), the lowest order 
equation ( ppj ) can be written in the form 

(^"^^ ^ dfiipSa/SUu) = Y 9fi!uMiaf3jf,u (46) 
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Substituting into equation (^Sf), this yields the resuh 

= '^'tiCiaCilS (^^^^^ - d^{pUy) (47) 

where the unit trace property of S^p has been used. In order to obtain an 
expression for Il^^j we use equations (^1]) and (|46|) to give, to first order in 
the velocity, 

in n, \ 

(48) 



0^') - -T 



where we have made the approximation, correct to 0(5^), that g'^^f' — gf^p- 
The term Srriia does not contribute to Il^^j, as can be seen from equation 



since it is included in the rest mass. It can further be shown that the 
contribution from the term Spiap is essentially zero. We therefore find 



c 



2 



(49) 



where we have substituted for the anisotropic relaxation parameter, r^. As- 
suming that the velocity tensors 

^c^i...a„ ~ ^^^jCjai Cia„ (50) 

i 

are isotropic up to sixth order, we recover the required form of the tensor 
coupling between the order tensor and the velocity gradient tensor. 

Using these results it can be shown that, apart from the term associ- 
ated with j3i, the required momentum evolution equation (||) is recovered 
when equation (|36| ) is used as the momentum forcing term. The term in 
/?! could be recovered by including a fourth order velocity product in the 
anisotropic relaxation parameter but this would require the velocity tensors 
to be isotropic to eighth order. However, for simplicity, the term in /3i is 
omitted in the results presented in this paper. 

The order tensor evolution equation is recovered from the two equations 
(|40|) and (^) by summing over the index i and adding the two resulting 
equations to give 

dtflipQafs) + dt^lipQafs) + d^{pQapU^) = ^ Xia/3 (51) 
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Equation (^Tj) combined with the forcing term gives the required evolu- 
tion equation (^). However, in order to obtain the result (|5TD it is necessary 
to suppress a contribution of the form 

(52) 

i 

at second order in the Chapman Enskog analysis and this is achieved using 
the term associated with the factor —2 in equation (0). 

The term arises because conservation of mass and momentum are 
alone insufficient to constrain all the first order elements of the distribution 
function. This arises because in order to recover the required macroscopic 
equations, the form of the equilibrium distribution function must satisfy 
equations ( ^61) . However, the limitations on the definitions of the moments 
of the distribution function lead to the result 

Ci-yQiap 7^ pSapu^ (53) 

i 

and as a consequence the argument of the derivative in (^) is non-zero. 



4.2 Choice of LB Parameters 



The following correspondence is found between the parameters of the LB 
algorithm and the parameters of the target scheme 
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/54 = ^-pci{l-2r,C) 

n P2 . 2 

P5 = -y + W^oC, 

= y + wroc^ (54) 

where the parameter ( = l + ?7o(l + '^/4) with d being the dimensionality 
of the LB scheme which is taken to be 2 for the results presented in Section 
d^). The viscosity set, {P^, P^, fii, fj,2}, of the Qian scheme is recovered from 
the parameter set, {tq, r/o, ''?2, /^i, /^2}, of the LB algorithm. We note from 
Qian( |2^) that Pq — P5 = fi2 and this is seen to be consistent with the last 



two of equations (|5^) . There is one free parameter in the LB scheme since we 
require to recover only four Qian parameters from five LB parameters. The 
free parameter is taken to be rjQ which is adjusted to place the LB scheme in 
a stable region of its parameter space, as is explained below. 
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The equations ( |5^ ) can be inverted to give 



4/iiC(2/35 + /i2) 

V2 = 



- ifiiiPi + pel) 



,2 



ro = (55) 

In the absence of any anisotropic terms, the last expression in equation (|55|) 
becomes p = c^(l — 2ro), which is equivalent to the standard result for 
the kinematic viscosity in an LBGK scheme. The relationship between the 
Qian parameters {/?4, /?5, /^i, /i2} and the standard Leslie Coefficients is given 
in Section |. 

We now establish an approximate criterion for the stability of the LB al- 
gorithm by considering a system which is evolving towards a uniform density 
and velocity distribution at all points. The algorithm involves repeated ap- 
plication of the collision operator, Mi^pj^u = M.i ^i-nd ignoring forcing terms, 
the non-equilibrium part of the distribution function is transformed at each 
time step according to the symbolic equation 

g(neq) ^ + ^) . ^("^9) (56) 

In order for the non-equilibrium part to decay to zero under successive ap- 
plications of equation (|56|), it is necessary that 

-l<{l + M)iaPj^u<l (57) 

for all choices of indices. This equation may be recast as the stability criteria 

- 2 < M^oiag) < \M^OffD^ag)\ < 1 (58) 



However, it should be remembered that the matrix M is a function of 
position, as a consequence of its dependence on the order tensor through the 
anisotropic relaxation parameter, Tj. Hence, in a given simulation, the condi- 
tions ( pSD cannot be guaranteed under all flow conditions. The criteria have 
been achieved in the implementation of the algorithm described in Section 
^ by explicit evaluation of all the non-zero elements of the collision matrix 
prior to running the algorithm and after making appropriate assumptions 
about the ordering in the system. The parameter 770 is then chosen to ensure 
that the conditions (|58[) are satisfied for all possible values of the director 
and scalar order parameter. 
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5 Analytical results from the Qian formalism 



In this section we quote without detailed proof some analytical results which 
can be derived from the Qian formalism and which will be used in Section |^ to 
validate the algorithm. Further discussion of the results and their significance 
is being prepared for a future publication. The results quoted in Section ^ 
have been derived for a system in which the director is confined to lie in a 
two dimensional plane and the order tensor for such a system may be written 
in the form 

Qaf3 = S{2nanf3 - (59) 

where the time and position dependent quantities S and tIq, are the order 
parameter and director. If the form (|59| ) is substituted into the free energy 
expression (P) and all the time and spatial derivatives are removed, we find 
the equilibrium Landau-deGennes free energy in the presence of a uniform 
magnetic field is given by 

Filc^ = + 475^ - ISnXaHl (60) 

where it is assumed that the director is parallel to the magnetic field. The 
equilibrium order parameter, which minimises equation (0), has the form 

Sh = S+ + (61) 



with 



1/3 



(62) 



where the free energy parameter, a, will be negative in the nematic phase 
and we choose the real solutions of Sq. In the limit the if — >■ 0, we have 



which is the equilibrium order parameter which effectively minimises the 
Landau-deGennes free energy in the absence of director gradients. Using 
equation ( |5P| ) with S set to Sq we may follow the arguments given in the 
Appendix (B) of Qian to determine the relationship between the Leslie coef- 
ficients of a material and the viscosity coefficients used in the Qian scheme. 
It is found for a system with the order tensor given by (p9|) 

Pi = ai/{ASi) P4 = (1/2) (2^4 + as + ae) 

(3, = a,/{2So) = a,/{2So) (64) 

Ail = {as - 02) / (SS^) H2 = {a2 + as) / {2So) 



14 



where 5*0 is the order parameter defined by equation (0). In simple shear 
fiow between two parallel plates and in the absence of any external field, 
equation for the evolution of the order tensor can be solved to find the 
steady state value of the angle of orientation and the order parameter. In the 
centre of the fiow, remote from the walls, the gradients in the order tensor 
may be ignored and the director is found to lie at an angle 6 with respect to 
the fiow velocity where 

cos(2^) = (65) 

and this is directly equivalent to the result found for the ELP equations, 
cos(26') = —71/72 ||10|- The order parameter is also modified in shear fiow 
and is given by the solution of the cubic equation 

ax^ + bx'^ + cx + d = Q (66) 

where x = S"^ and 

a = 102472 b = 256^7 , . 

c = 16(a2 + {iJ,idyU^)^) d = -{^i2dyU.^f 

This equation may be solved and to first order in the velocity gradient the 
order parameter is found to be 

b = bQ + - . Idyti^l (68) 

o\Olp\ 

and hence the order parameter increases as a consequence of fiow alignment. 
Finally we consider the Miesowicz viscosities in a system with a variable 
order parameter. The stress tensor evolution equation can be solved for a 
simple shear fiow in the presence of a strongly aligning external field. If the 
director is given by {ux, Uy], the associated Miesowicz viscosity, ?7e//, is given 

by 

Veff = {S'o[2a4 + as + tte] + 5'2[a3 - "2 + 4ain^.nJ] 

+SSo[aQ{2nl - 1) + (^2 + «3)(n^ - nj) + a5(2nj - 1)]} (69) 

where S is the order parameter in the fiow and So is the equilibrium value 
of the order parameter defined in equation (|63|). In the limit that the order 
parameter is fixed {ie S = Sq), equation ([H5|) reduces to the standard results 



Vpar = ^(a4 + a3 + «6) 
Vperp = ^{0!4-a2 + a5) (70) 
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where rjpar and r/perp correspond to the director being ahgned parallel, n = 
{1,0} and perpendicular, n = {0, 1}, to the flow. In the presence of both a 
velocity gradient and a strong aligning field the change in the order parameter 
is dominated by the aligning field and is hence given to a good approximation 
by equations (|6T|) and (B^). 



6 Results 



A number of simulations were undertaken in order to validate the method 
described in the previous section. All the results were obtained using a 
D2Q13 lattice. 

In the absence of flow and in the presence of periodic boundary conditions, 
the order parameter is found to be consistent with the value predicted by 
equation ( |63D to machine accuracy. If an external magnetic field is applied 
through the introduction of a molecular field term of the form (|T^) , it is 
found that the director becomes fully aligned with the external field and that 
the order parameter follows equation ( |6lD with an accuracy of better than 
0.1%. These two results confirm that the free energy correctly controls the 
dynamical equations through the associated molecular field. In principal, 
therefore, the effects of temperature could be introduced into the model 
through the coefficients in the Landau-deGennes free energy. 

The angle of alignment of the director in a shear flow in the absence of 



an external magnetic field is given by equation (|65D. In order to test that the 
technique recovered this result, simulations were run with the free energy 
parameters, {ap = — 0.0512, /3f = 0,7i? = 0.01), chosen to give an equi- 
librium order parameter of 0.8 and the parameters, (Li = 0.001,^2 = 0), 
selected to give nematic elastic behaviour equivalent to the one constant 
approximation. The LB parameters were chosen to be, (tq = 1.1000, ?7o = 
0.3036, r/2 = 0.1565, /ii = 0.3823, s/ia = -1.261, p = 1.8). If the order pa- 
rameter 5*0 is assumed to be 0.8, these values recover the viscosity ratios, 
(02/04 = —0.9556,03/04 = —0.0144,05/04 = 0.5565,05/04 = —0.4135), 



which are consistent with MBBA at 25° C |E3 



With this choice of parameters, the angle of rotation of the director is 
predicted to be 7° and this is recovered to an accuracy of less than 1%. The 
parameters /ii and fj,2 were then adjusted to give a range of values of angle of 
alignment and the results are shown in Figure (|Tp together with the expected 
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curve, equation (0), which is shown as a continuous curve with the value 
of the order parameter taken from the simulation data. The change in the 
order parameter predicted by equation ( pE| ) is also recovered to within less 
than 0.3%. 

As a final test of the algorithm, the Miesowicz viscosities were measured 
for simulations in which the director orientation was controlled by a strong 
magnetic field. Using the parameters appropriate for MBBA the expected 
ratio of the viscosity with the director aligned parallel and perpendicular to 
the direction of flow is 5.03 and this value was recovered with an error of 
0.1%. The expected form of the viscosity is given by equation ( p^ and it 
can be shown that this implies that the viscosity should be a linear function 
of I/tq. Figure shows that the expected linearity is observed within the 
simulations; the slope to intercept ratio for each curve is in error by less than 
0.2% 



7 Conclusions 



In this paper we have presented a generalised lattice Boltzmann scheme which 
recovers the tensor order parameter equations for a nematic liquid crystal 
proposed by Qian and Sheng The generalised method is based on a 



single tensor density from which all the macroscopic quantities can be recov- 
ered. A Chapman Enskog analysis demonstrates that the algorithm recovers 
the target macroscopic equations and test simulations demonstrate that the 
method correctly recovers the evolution of the director, the order parameter 
and the velocity gradients in the presence of shear and magnetic fields. 

The method will generalise straightforwardly to three dimensions and 
work is currently in progress to extend the scheme to model a mixture of a 
isotropic and nematic fluids. 
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Figure 2: Viscosity in shear flow as a function of I/tq. The upper (lower) 
line is for director perpendicular (parallel) to the flow direction. 
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